% Taking the firm-side calibration as given, this script does the following:
% 1) It loads the data corresponding to the agent's behaviour under an
% alternative specification of the withdrawal penalty schedule. 
% For example: The current script loads the data where the penalty is equal
% to 0.25*the baseline penalty. 
% 2) Given the firm-side calibration, it calculates the equilibrium profit
% level at the laissez-faire rate of return and then identifies the new
% equilibrium by looking for a contract on the modified iso-profit which
% maximises the agent's perceived utility. Doulbe-check of no profitable
% price-based deviations is built in.
% 3) It compares the outcomes to laissez-faire. 


clear
load('recalibrated_data_zpen025_with_eval.mat')

fileID = fopen('Policy_3_penalty_2021.txt','wt')



rho = mydata.rho;
delta = mydata.delta;
effhh_size_ = mydata.effhh_size;
avgIncome = mydata.avgIncome;
meanhhs = mydata.meanhhs;
R = mydata.R;
zpen_eval = mydata.zpen_eval;
alpha = mydata.alpha;
survival = mydata.survival;
age_ = mydata.age;

THETAtarget = 0.05;
totalcosttarget = 0.005;
adminratiotarget = 0.50;
markuptarget = 0.2;

suffix = strcat('avgZ',num2str(THETAtarget*10000),'_naif');
avgZ_naif_target = mydata.(suffix);
suffix = strcat('fmax_',num2str(THETAtarget*10000));
fmax_target = mydata.(suffix);


% Set the calibrated values:
gamma1 = 0.5;
c1 = adminratiotarget*totalcosttarget*mean(avgZ_naif_target) / mean(avgZ_naif_target.^(gamma1));
gamma2 = 5.75;
c2 = (1-adminratiotarget)*totalcosttarget*mean(avgZ_naif_target) / (1+THETAtarget-R)^(gamma2);
%phi = 0.35;
xi = 0.481482;

% And for later comparisons with laissez-faire:
lf_profit = 15211.219049;
lf_fee = 1919.75;
lf_actual_utility = 654.673018;
lf_savings = 501488;


% 1) Given the Hotelling parameter xi, what fee level is not competed away
% at Rz = 5% ?

fee_found = 0;

if rho == 1
    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
else
    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
end

suffix = strcat('avgZ',num2str(THETAtarget*10000),'_TC');
avgZ_TC = mydata.(suffix);
suffix = strcat('avgX',num2str(THETAtarget*10000),'_TC');
avgX_TC = mydata.(suffix);
suffix = strcat('avgTC',num2str(THETAtarget*10000),'_TC');
avgTC_TC = mydata.(suffix);

bequest = zeros(1,71);
for t = 1:71
    cons_bequest = (R-1)*(max(0,avgX_TC(t)+(1-zpen_eval(t))*avgZ_TC(t))) + avgIncome;
    if rho == 1
        yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    else
        yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    end
    bequest(t) = alpha*(utility_withbeq - utility_nobeq);
end

fee_at_target = lf_fee;
while fee_found == 0
    % 1a) Calculate utility from the (supposed) equilibrium contract
    avgTC_fee = avgTC_TC - fee_at_target*ones(1,71);
    
    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    lifetime_equilibrium = instant(1);
    for t = 2:71
        lifetime_equilibrium = lifetime_equilibrium + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
    
    gridlength = 11;
    gridjump = 0.05;

    ProfitPVvec = zeros(1,gridlength);
    expected_profits = zeros(1,gridlength);
    
    for j = 1:gridlength
    
        % First, profit from an interaction
        fee = (1 - floor(gridlength/2)*gridjump + (j-1)*gridjump)*fee_at_target;
    
        THETA = THETAtarget;
        avgZ_ = avgZ_naif_target;
    
        costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);

        profits = ones(1,71)*fee - costs;

        ProfitPV = profits(1);
        for t = 2:71
            ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
        end
    
        ProfitPVvec(j) = ProfitPV;
    
        % Second, probability of interaction
        avgTC_fee = avgTC_TC - fee*ones(1,71);
    
        instant = zeros(1,71);
        for t = 1:71
            if rho == 1
                instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
            else
                instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
            end
        end
    
        lifetime_fee = instant(1);
        for t = 2:71
            lifetime_fee = lifetime_fee + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
        end
    
        probabilityraw = 0.5 + (lifetime_fee - lifetime_equilibrium)/(2*xi);
        probability = max(0,min(1,probabilityraw));
    
        expected_profits(j) = ProfitPV*probability;
    end
    
    if max(expected_profits) == expected_profits(ceil(gridlength/2))
        fee_found = 1;
        fprintf(fileID, ['************** \r\n']);
        fprintf(fileID, ['For xi = %f, the calibrated fee of %f is the new Hotelling equilibrium \r\n'],[xi,fee_at_target]);
    else
        [maximised_profits, id_max] = max(expected_profits);
        fee_at_target = (1 - floor(gridlength/2)*gridjump + (id_max-1)*gridjump)*fee_at_target;
    end
end


% 2) Calculate new iso-fees.

% 1c) Profit level at the target
THETA = THETAtarget;
avgZ_ = avgZ_naif_target;
    
costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);

profits = ones(1,71)*fee_at_target - costs;

maxprofit = profits(1);
for t = 2:71
    maxprofit = maxprofit + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
end


% 1d) Need to re-calculate the iso-fees
for THETA = mydata.THETA_grid
    if abs(THETA - THETAtarget) < eps
        suffix = strcat('isofee_',num2str(THETA*10000));
        workspace.(suffix) = fee_at_target;
    else
        suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
        avgZ_ = mydata.(suffix);
                        
        costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);
                        
        fee = 100;
        profits = ones(1,71)*fee - costs;
        ProfitPV = profits(1);
        for t = 2:71
            ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
        end
                    
        while maxprofit > ProfitPV
            fee = fee + 6;
            profits = ones(1,71)*fee - costs;
            ProfitPV = profits(1);
            for t = 2:71
                ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
            end
        end
                        
        suffix = strcat('fmax_',num2str(THETA*10000));
        fmax = mydata.(suffix);
                        
        suffix = strcat('isofee_',num2str(THETA*10000));
        workspace.(suffix) = min(fee-3,fmax);
    end
end

                    
% 1e) Calculate profits for each Rz (just to double-check the iso-profit condition)
fprintf(fileID, ['**************\r\n']);
for THETA = mydata.THETA_grid
    suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
    avgZ_ = mydata.(suffix);
    
    costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);
    
    suffix = strcat('isofee_',num2str(THETA*10000));
    fee = workspace.(suffix);
         
    profits = ones(1,71)*fee - costs;
    ProfitPV = profits(1);
    for t = 2:71
        ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
    end
             
    fprintf(fileID, ['For Rz = %f, the PV of profits (on the iso-profit curve) is %f \r\n'],[THETA,ProfitPV]);
end


% 2) Given the new markup at Rz = 5% and recalculated iso-fees, which
% contract maximises perceived utility and is therefore offered in
% equilibrium?
fprintf(fileID, ['**************\r\n']);

if rho == 1
    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
else
    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
end
                
perceived_lifetime_iso_grid = [];
for THETA = mydata.THETA_grid
    suffix = strcat('avgZ',num2str(THETA*10000),'_TC');
    avgZ_TC = mydata.(suffix);
    suffix = strcat('avgX',num2str(THETA*10000),'_TC');
    avgX_TC = mydata.(suffix);
    suffix = strcat('avgTC',num2str(THETA*10000),'_TC');
    avgTC_TC = mydata.(suffix);
                    
    suffix = strcat('isofee_',num2str(THETA*10000));
    isofee = workspace.(suffix);
    avgTC_TC = avgTC_TC - ones(1,71)*isofee;

    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_TC(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_TC(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    bequest = zeros(1,71);
    for t = 1:71
        cons_bequest = (R-1)*(max(0,avgX_TC(t)+(1-zpen_eval(t))*avgZ_TC(t))) + avgIncome;
        if rho == 1
            yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        else
            yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        end
        bequest(t) = alpha*(utility_withbeq - utility_nobeq);
    end
    
    lifetime = instant(1);
    for t = 2:71
        lifetime = lifetime + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
                    
    perceived_lifetime_iso_grid = [perceived_lifetime_iso_grid lifetime];
    
    fprintf(fileID, ['For Rz = %f and the associated iso-fee, the PERCEIVED lifetime utility is %f \r\n'],[THETA,lifetime]);
                    
    if abs(THETA - THETAtarget) < eps
        lifetime_at_target = lifetime;
    end
end

[market_benchmark, id_market] = max(perceived_lifetime_iso_grid);
fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['Now, Rz =  %f maximises perceived utility and is offered in market equilibrium \r\n'],[mydata.THETA_grid(id_market)]);



% If Rz offered in new equilibrium differs from laissez-faire, need to
% double-check no profitable price deviations condition:
% Verify that there is no incentive to price-deviate from
% iso_fee(id_market) by +/- 25%: 
THETA = mydata.THETA_grid(id_market);

suffix = strcat('isofee_',num2str(THETA*10000));
fee_supp_eq = workspace.(suffix);
suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
avgZ_ = mydata.(suffix);

suffix = strcat('avgZ',num2str(THETA*10000),'_TC');
avgZ_TC = mydata.(suffix);
suffix = strcat('avgX',num2str(THETA*10000),'_TC');
avgX_TC = mydata.(suffix);
suffix = strcat('avgTC',num2str(THETA*10000),'_TC');
avgTC_TC = mydata.(suffix);

if rho == 1
    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
else
    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
end

bequest = zeros(1,71);
for t = 1:71
    cons_bequest = (R-1)*(max(0,avgX_TC(t)+(1-zpen_eval(t))*avgZ_TC(t))) + avgIncome;
    if rho == 1
        yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    else
        yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    end
    bequest(t) = alpha*(utility_withbeq - utility_nobeq);
end

% Calculate perceived utility at current eqm candidate
avgTC_fee = avgTC_TC - fee_supp_eq*ones(1,71);    
instant = zeros(1,71);
for t = 1:71
    if rho == 1
        instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
    else
        instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
    end
end
lifetime_new_eq = instant(1);
for t = 2:71
    lifetime_new_eq = lifetime_new_eq + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
end

% Verify no profitable price deviation: 
gridlength = 11;
gridjump = 0.05;

ProfitPVvec = zeros(1,gridlength);
expected_profits = zeros(1,gridlength);

for j = 1:gridlength
    % First, profit from an interaction
    fee = (1 - floor(gridlength/2)*gridjump + (j-1)*gridjump)*fee_supp_eq;
    costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);
    profits = ones(1,71)*fee - costs;
    ProfitPV = profits(1);
    for t = 2:71
        ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
    end
    ProfitPVvec(j) = ProfitPV;
    
    % Second, probability of interaction
    avgTC_fee = avgTC_TC - fee*ones(1,71);
    
    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    lifetime_fee = instant(1);
    for t = 2:71
        lifetime_fee = lifetime_fee + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
    
    probabilityraw = 0.5 + (lifetime_fee - lifetime_new_eq)/(2*xi);
    probability = max(0,min(1,probabilityraw));
    
    expected_profits(j) = ProfitPV*probability;
end

[max_profit,adjustment_id] = max(expected_profits);
optimal_adjustment = (1 - floor(gridlength/2)*gridjump + (adjustment_id-1)*gridjump);
fee_eq = optimal_adjustment*fee_supp_eq;
fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['For xi = %f and Rz = %f, there is no profitable price-based deviation for fee = %f  \r\n '],[xi,THETA,fee_eq]);
fprintf(fileID, ['(Adjustment (if any) with optimal multiplier of  %f implemented)  \r\n '],[optimal_adjustment]);

% Note: Need to multiply max_profit by 1/0.5 to correct for eqm
% probability of interaction:
max_profit = max_profit*2;


% 1f) Calculate which of the contracts lying on the iso-profit maximises
% ACTUAL consumer welfare, to find Pareto-optimum. 
% Note: A different Rz may be Pareto-efficient when the withdrawal penalty
% schedule changes!
fprintf(fileID, ['**************\r\n']);
if rho == 1
    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
else
    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
end

actual_lifetime_iso_grid = [];
for THETA = mydata.THETA_grid
    suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
    avgZ_ = mydata.(suffix);
    suffix = strcat('avgX',num2str(THETA*10000),'_naif');
    avgX_ = mydata.(suffix);
    suffix = strcat('avgTC',num2str(THETA*10000),'_naif');
    avgTC_ = mydata.(suffix);
                    
                    
    suffix = strcat('isofee_',num2str(THETA*10000));
    isofee = workspace.(suffix);
    avgTC_ = avgTC_ - ones(1,71)*isofee;

    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    bequest = zeros(1,71);
    for t = 1:71
        cons_bequest = (R-1)*(max(0,avgX_(t)+(1-zpen_eval(t))*avgZ_(t))) + avgIncome;
        if rho == 1
            yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        else
            yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        end
        bequest(t) = alpha*(utility_withbeq - utility_nobeq);
    end
    
    lifetime = instant(1);
    for t = 2:71
        lifetime = lifetime + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
                    
    actual_lifetime_iso_grid = [actual_lifetime_iso_grid lifetime];
    %fprintf(fileID, ['For Rz = %f and the associated iso-fee, the ACTUAL lifetime utility is %f \r\n'],[THETA,lifetime]);
                    
end

[fb_utility, id_fb] = max(actual_lifetime_iso_grid);
fprintf(fileID, ['The efficient contract has Rz =  %f \r\n'],[mydata.THETA_grid(id_fb)]);


% Impact on firms profits
profitchange = (maxprofit - lf_profit)/lf_profit;

fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['Relative to laissez faire, the change in firms profits is %f percent\r\n'],[100*profitchange]);


% 2b) Calculate the ACTUAL consumer welfare assocviated with the new
% equilibrium contract 

fprintf(fileID, ['**************\r\n']);
if rho == 1
    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
else
    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
end

actual_lifetime_iso_grid = [];
for THETA = mydata.THETA_grid(id_market)
    suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
    avgZ_ = mydata.(suffix);
    actual_savings_at_retirement = avgZ_(45);
    suffix = strcat('avgX',num2str(THETA*10000),'_naif');
    avgX_ = mydata.(suffix);
    suffix = strcat('avgTC',num2str(THETA*10000),'_naif');
    avgTC_ = mydata.(suffix);
    
    suffix = strcat('isofee_',num2str(THETA*10000));
    fee = workspace.(suffix);
                    
    avgTC_ = avgTC_ - ones(1,71)*fee;

    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    bequest = zeros(1,71);
    for t = 1:71
        cons_bequest = (R-1)*(max(0,avgX_(t)+(1-zpen_eval(t))*avgZ_(t))) + avgIncome;
        if rho == 1
            yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        else
            yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        end
        bequest(t) = alpha*(utility_withbeq - utility_nobeq);
    end
    
    lifetime = instant(1);
    for t = 2:71
        lifetime = lifetime + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
                    
    %actual_lifetime_iso_grid = [actual_lifetime_iso_grid lifetime];
    fprintf(fileID, ['For Rz = %f and the associated iso-fee, the ACTUAL lifetime utility is %f \r\n'],[THETA,lifetime]);
                    
end


% What is the CE welfare loss?
utility_multiplied = lifetime;
efficient = lf_actual_utility;
multiplier = 1.00;

if utility_multiplied < efficient
    while utility_multiplied < efficient 
        instant_multiplied = zeros(1,71);
        for t = 1:71
            if rho == 1
                instant_multiplied(t) = effhh_size_(t) * log((multiplier*avgTC_(t))/effhh_size_(t));
            else
                instant_multiplied(t) = effhh_size_(t) * (((multiplier*avgTC_(t))/effhh_size_(t))^(1-rho)-1) / (1-rho);
            end
        end
    
        utility_multiplied = instant_multiplied(1);
        for t = 2:71
            utility_multiplied = utility_multiplied + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant_multiplied(t) + (1-survival(t))*bequest(t));
        end
        multiplier = multiplier + 0.0001;
    end

    fprintf(fileID, ['The CE welfare LOSS relative to laissez faire is %f percent\r\n'],[100*(multiplier-1.0001)]);

else
    while utility_multiplied > efficient
        instant_multiplied = zeros(1,71);
        for t = 1:71
            if rho == 1
                instant_multiplied(t) = effhh_size_(t) * log((multiplier*avgTC_(t))/effhh_size_(t));
            else
                instant_multiplied(t) = effhh_size_(t) * (((multiplier*avgTC_(t))/effhh_size_(t))^(1-rho)-1) / (1-rho);
            end
        end
        
        utility_multiplied = instant_multiplied(1);
        for t = 2:71
            utility_multiplied = utility_multiplied + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant_multiplied(t) + (1-survival(t))*bequest(t));
        end
        
        multiplier = multiplier - 0.0001;
        
    end
    fprintf(fileID, ['The CE welfare GAIN relative to laissez faire is %f percent\r\n'],[100*(1-0.0001-multiplier)]);
end

        
% Finally, what is the impact on retirement savings?
savings_impact = (actual_savings_at_retirement - lf_savings)/ lf_savings;
fprintf(fileID, ['Relative to laisseiz faire, the savings at retirement change by %f percent \r\n'],[100*(savings_impact)]);    



    
    
                    
                    
